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Correlation effects in solids from first principles 
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1-1-4 Higashi, Tsukuba, Ibaraki 305-0046, Japan 



Abstract 

First principles calculations of bandstructures of crystals are usually based on one-particle 
theories where the electrons are assumed to move in some effective potential. The most 
commonly used method is based on density functional theory within the local density ap- 
proximation (LDA). There is, however, no clear justification for interpreting the one-particle 
eigenvalues as the bandstructure. Indeed, the LDA failure to reproduce the experimental 
bandstructure is not uncommon. The most famous example is the bandgap problem in 
semiconductors and insulators where the LDA generally underestimates the gaps. 

A rigorous approach for calculating bandstructures or quasiparticle energies is provided 
by the Green function method. The main ingredient is the self-energy operator which acts like 
an effective potential but unlike in the LDA, it is nonlocal and energy dependent. The self- 
energy contains the effects of exchange and correlations. An approximation to the self-energy 
which has proven fruitful in a wide range of materials is the so-called GW approximation 
(GWA). This approximation has successfully cured the LDA problems and has produced 
bandstructures with a rather high accuracy. For example, bandgaps in s-p semiconductors 
and insulators can be obtained typically to within 0.1-0.2 eV of the experimental values. 

Despite its success, the GWA has some problems. One of the most serious problems is its 
inadequacy to describe satellite structures in photoemission spectra. For example, multiple 
plasmon satellites observed in alkalis cannot be obtained by the GWA. Recently, a theory 
based on the cumulant expansion was proposed and shown to remedy this problem. Apart 
from plasmon satellites which are due to long-range correlations, there are also satellite 
structures arising from short-range correlations. This type of satellite cannot be described 
by the cumulant expansion. A t-matrix approach was proposed to account for this. 

Although traditionally the Green function method is used to calculate excitation spectra, 
groundstate energies can also be obtained from the Green function. Recent works on the 
electron gas have shown promising results and some approaches for calculating total energies 
will be discussed. 
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1 Introduction 

The electronic structures of molecules and solids are now routinely calculated using the Kohn- 
Sham (KS) density functional theory (DFT) [52, 61] within the local density approximation 
(LDA). In this theory, one solves a single-particle Hamiltonian with an effective local potential 
containing the Hartree potential Vh and an exchange-correlation potential V xc : 



The resulting one-particle eigenvalues e\ are often interpreted as one-particle excitation energies 
or in the case of a crystal, as quasiparticle energies (bandstructures) measured in photoemission 
experiments. There is no clear theoretical justification for this, except for the highest occupied 
states. For sp systems, it is found that the bandstructure is often in rather good agreement 
with experiment. It is remarkable that the LDA, being a simple approximation, can go a long 
way in accounting for the electronic structures of many materials, However, there are serious 
discrepancies. The most well-known of these is the band-gap problem in semiconductors (Si, 
GaAs, Ge, etc.) where the LDA bandgpas are systematically underestimated. This is illustrated 
in figure 1. The discrepancies become worse in strongly correlated 3d and 4/ systems. In the 
so called Mott-Hubbard insulators of transition metal oxides, the LDA band gap is much too 
small. In some cases, the LDA even gives qualitatively wrong results. For example, the Mott- 
Hubbard insulator CoO and the undoped parent compound of the high T c material La2CuO"4 are 
predicted to be metals whereas experimentally they are insulators. Apart from these problems, 
LDA sometimes overestimates the valence bandwidth, for example in Na and Ni. 

A proper way of calculating single-particle excitation energies or quasiparticle energies [64, 
65] is provided by the Green function theory [38, 39]. It can be shown that the quasiparticle 
energies Ei can be obtained from the quasiparticle equation: 



The non-local and energy dependent potential S, or the self-energy, contains the effects of 
exchange and correlations. It is in general complex with the imaginary part describing the 
damping of the quasiparticle. It is then clear that the different single-particle theories amount 
to approximating the self-energy operator S. Thus we can think of the V xc in DFT as a local 
and energy independent approximation to the self-energy that gives the correct ground state 
density. Approximating S by the exchange operator 



results in the Hartree-Fock approximation. 

Improving the LDA means finding corrections to the self-energy beyond V r I ^ DA . There are a 
number of attempts, amomg these are Generalized Gradient Approximations (GGA) [67, 16, 17, 
18, 101, 95, 88], Self-Interaction Correction (SIC) [26, 71, 108, 87, 100, 102, 7], LDA+U [4, 5, 
68, 93, 94, 6], Optimized Effective Potential (OEP) Method (Exact Exchange + Correlations) 
[103, 62, 63, 22, 23], and more recently Dynamical Mean-Field theory (DMFT) (for a review 
and references, see [40]). This note describes the GW approximation (GWA) [46] which has 
been found to be successful in accounting for quasiparticle energies for a wide range of systems 
from atoms to solids. The GWA may be thought of as a generalization of the Hartree-Fock 
approximation (HFA) but with a dynamically screeneed interaction. Recent development in 
going beyond the GWA is also presented. These include the cumulant expansion method and 
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the t-matrix approach. In addition, applications of the Green function theory to total energy 
calculations and the issue of self-consistency will be discussed. 

This note also describes practical implementations of the GWA. Conventional methods for 
calculating the self-energy use plane waves as basis functions which are suitable for sp systems. 
However, applications to 3d and 4f systems using plane-wave basis are unsuitable due to a large 
number of plane waves needed to describe the localized 3d or 4f states. A method for treating 
systems with localized states is described. The method is based on the linear muffin-tin orbital 
(LMTO) approach [3] and basis functions needed to describe the screened Coulomb interaction 
and the self-energy are products of the LMTO's. 

2 Theory 

2.1 Photoemission experiments measure quasiparticle energies 

In a photoemission experiment, photons are used to excite electrons out of a crystal leaving 
holes behind and giving information about the occupied states. In an inverse photoemission 
experiment (BIS), electrons are sent into the crystal to probe unoccupied states. Photoemission 
processes are usually interpreted in terms of the three-step model [34, 81, 89, 20] consisting of 

• optical excitation of an electron; 

• its transport through the solid with the possibility of inelastic scattering with other elec- 
trons; 

• the escape through the surface into the vacuum. 

To illustrate the basic principle, here we consider a simple case of normal emission where photo- 
electrons are emitted normal to the sample surface. The more general case is considerably more 
complicated and for a more detailed discussion, we refer to Ref. [92]. 
Using conservation of energy 

E (N)+l; = E f {N - l,k) + A; 2 /2 (4) 

where Eo(N) is the groundstate energy of the solid, to is the photon energy, E/(N — l,k) is an 
excited state of the (N-l) electrons, and k 2 /2 is the kinetic energy of the outgoing photoelectron. 
For a given crystal wave vector k , the photoelectron spectrum usually shows peaks at well 
defined energies w - k 2 /2 = E/(N - 1, k) - Eq(N) which may be interpreted as the quasiparticle 
excitations and they can be shown to correspond to peaks in the imaginary part of the Green 
function. When the positions of the peaks are plotted as a function of k , we obtain what is 
called the bandstructure. 

^From the Fermi Golden Rule, the probability per unit time of emitting an electron of 
momentum k is 

I (k, w) = 2tt £ \{N - 1, /, k| A ■ p\N) \ 2 S(u - k 2 /2 - e f ) (5) 
/ 

where A is the electromagnetic vector field, £/ — Ef(N — 1) - Eo(N). Since the photoelectron 
usually has a high energy we can approximate ("sudden approximation") 

|7V-l,/,k>«4|7V-l,/), c k |iV)«0 (6) 

In second quantized form A • p = J2im ^im c ] c m- Using c^cj + cjc^ = 6^ we then have 

/(k, U ) = 2n ]T \{N -1J\J2 ^ m c m \N)\ 2 S(uj - k 2 /2 - Ef) (7) 

f m 
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which can be rewritten 

I{k,u) = 2tt 2 A km A mm ,(w - k 2 /2)A km , (8) 
mm' 

where 

Amm'(u) = £{iV - l,/|c m |iV){iV|c m ,|iV - l,/)<5(c - €/) (9) 

/ 

The quantity A will be shown later to be proportional to the imaginary part of the Green 
function, assuming the matrix elements Ak m are constant. 

2.2 The Green function and the self-energy 

To describe the photoelectron spectra measured in photemission and inverse photoemission 
experiments we need a quantity that the propagation of a hole or an added electron. A suitable 
quantity is the one-particle Green function Refs. [36, 78, 59] which is defined as: 

iG(x,x') ■= {N\T[ij>{x)ft(x')}\N) 

f (N\ip{x) ft(x')\N) for t > t' (electron) , . 

\ -(N\ft(a>) j>(x)\N) for t < t' (hole) ( ' 

We have used a notation x = (r, t). T is the time ordering operator, and |iV) is the groundstate 
of TV electrons. Thus, for t > t', the Green function is the probability amplitude that an electron 
added at x' will propagate to x, and for t' > t, the probability amplitude that a hole created at 
x will propagate to x'. A possible spin flip can be incorporated in the definition of G. 
^From the Green function we can obtain: 

• The expectation value of any single-particle operator in the ground state. 

• The ground state energy 

• The excitation spectrum 

The field operators are defined in the Heisenberg representation and they satisfy the equation 
of motion 

i^(x) = $(x),A] (11) 

where the Hamiltonian is given by 

H = j d 3 r ft(x)h (r)$(x) 

+ \f d 3 rd 3 r' ij>\x)ft (x?)v(r - r')V>(x')^(x) (12) 

ho is the kinetic energy operator plus any single-particle operator such as an external field. 
Evaluating the above commutator, the equation of motion for the Green function is then 

d 1 

+ i f d 3 n v{v-v x ) (N\T[^\r u t)^(r u t)^(r,t)^Hr',t')]\N) 
= S(x-x') (13) 
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The quantity 

(N\T[^(r u t)^( ri ,t)^(r,t)^H^^)]\N) 
is a special case of the two-particle Green function defined by 

G 2 (l,2,3,4) = ■(i) 2 (N\T[fal)$(3)ft(4)ft(2)\N) (14) 

describing the propagation of two particles from 2, 4 to 1,3 where we use a short- hand notation 
1 = x\. The two-particle Green function, in turn, satisfies an equation of motion involving a 
three-particle Green function, forming a hierarchy of equations. To break the hierarchy, the 
two-particle Green function is expressed in terms of the one-particle Green functions. Formally, 
we introduce the mass operator M(x,x') such that 

jdx l M { x,x l )Gix u x<) = 

-ij d 3 n «(r-n) {N\T[^(r l ,t)^{r 1 ,t)^(r,t)^{T f ,lf)]\N) (15) 



so that 



G(x,x') - j dxiM{x,xi)G(x u x') = S(x - x') (16) 



As will be shown later, the mass operator contains the Hartree potential and the self-energy is 
defined to be the mass operator without the Hartree potential, thus containing the effects of 
exchange and correlation only: 

S = M - V H (17) 
Fourier transforming Eq. (16) and using the definition of the self-energy we get 

[w-Mr)-Vff(r)]G(r,r' ) w)-|dr 1 E(r,r 1 ,tj)G(r 1 ,r',a;)==*(r-r') (18) 

We define a noninteracting Green function Goasa solution to Eq. (18) with E = 0. Noting 
that ui - ho — Vh — Gq l , the exact Green function G satisfies the Dyson equation 



Go + Go EG 
1 



Go 1 - E 



(19) 



3 Spectral representation of the Green function 



The spectral representation of the Green fucntion is obtained from the definition in Eq. ( 10) 
by noting that 

i>{r,t) = exp{iHt) V>(r,0) exp(-iHt) (20) 

Inserting a complete set of N ± 1-particle eigenstates of H and performing the Fourier transform 
give 

r<r r' , A = V MOW + V PiW) m ) 

K ' ' ; ^ u, + Ei(N - I) - E {N) - i6 yu-E 3 (N + l) + Eo(N) + i6 



where 



hi(r) = <JV-l,#(r,0)|iV) . (22) 
Pj (r) = {N\^{r,0)\N + lJ) (23) 
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\N ± 1, i) is the ith eigenstate of the N ± 1 electrons with eigen energy Ei(N±l) and Eq(N±1) 
is the groundstate energy of the N ± 1 electrons. 

We can write the Green function in a more compact form 

G(r,r',w)=/ dJ , .' + drf /. (24 

7-oo tjJ — U> ' — 10 U) — U)' + 10 

The spectral function A (density of states) is proportional to the imaginary part of G 

A(r, r',u>) = -— Im G(r, r', w)sgn(ai - fi) (25) 

7T 

and it is given by 

4(r,r',w) - ^/iW/'fr'^-Zi + efiV-M)] 

i 

+ Y,9i(T)9?(r')6[cj-ti-e{N + l,i)] (26) 

X 

d(N ± 1) is the excitation energy of the N ± 1 electrons: 

d(N ± 1) = £i(AT ± 1) - £ (iV ± 1) (27) 

which is positive and the quantity fi is the chemical potential 

ft = £(JV + 1) - £(iV) = £(JV) - E(JV - 1) + 0(l/iV) (28) 

The poles of the Green function are therefore the exact excitation energies of the N±l electrons. 

As an example, for a noninteracting case the many-body states become single-Slater deter- 
minants. The quantities hi and pi are replaced by the one-particle wavefunctions satisfying 

Hotpi = £iipi (29) 

and the excitation energies are replaced by e,. The noninteracting Green function is then given 

Go(r , r>) = £«M<£l/!fW^ m 

i 3 

The energy is measured with respect to the chemical potential /i. 
3.1 Quasiparticles 

For atoms and molecules, the excitation energies are discrete and the poles of the Green function 
are consequently discrete. In crystals, the excitation energies become continuous and the poles 
of the Green function form a branch cut. It is then not meaningfull to speak about individual 
poles. It is here the concept of quasiparticle comes in and it may be understood as follows. 
^From Eq. ( 19), the spectral function A is schematically given by 

A(u>) = I^IImGiMI 

7T . 
I 

1^ |Im Sj (a;)) 

n ^ |w - £i - Re Ei{co)\ 2 + |Im £j(u;)| 2 { ' 

where Gj is the matrix element of G in an eigenstate ipi of the noninteracting system Hq. In 
a crystal, the state label i corresponds to the wavevector k and band index n. A is usually 
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peaked at each energy Ei = - Re £i(i?j) provided ImS(jE'j) is small. This peak is called a 
quasiparticle peak with a life-time given by |Im Ei(Ei)\. The life-time is a consequence of the 
fact that quasiparticle excitations are not eigenstates of the interacting Hamiltonian. 

The weight of the quasiparticle (strength of the Lorentzian) can be obtained by Taylor 
expanding ReS around the quasiparticle energy: 

ReSM = ReS(^) + (u - Ei) d **f(Ei) + 

ou 



We obtain for u close to Ei 



|w - £;| 2 + |Im Ei(o;)| 2 



where 

1 dRe-Zi(Eiy- ] 

— 



< 1 (34) 



^From the definition of G we have 

j du>A(r, r'; u) = (5(r - r') (35) 



or 



J du)Ai{u) = 1 (36) 



Since Z{ is usually less than one, this means the rest of the weight must be distributed at 
other energies. At some other energies u p , the denominator may be small and A(u) p ) could also 
show peaks or satellite (incoherent) structure which can be due to plasmon excitations or other 
collective phenomena. This is in contrast to the nonineracting case where the quasiparticle peak 
is just a delta function with an infinite life-time and without any incoherent features. 

^From the classical theory of the Green functions a solution to equation (18) can be written 

as 

C (ry,^E *' ( ^>!-/ r >> (37) 

where ^ are solutions to 

# (r)*i(r,uO + j d*rV{T,T\u)<!! i {r',u>)=E i {u)y i {T,L>) (38) 

We define a quasiparticle wave function with energy Ei as a solution to 

#o(r)* 2 (r) + j dr'S(r,r',^)^(r') - £^(r) (39) 

The eigenvalues Ei are in general complex and the quasiparticle wavefunctions are not in general 
orthogonal because S is not Hermitian but both the real and imaginary part of E are symmetric. 

3.2 Derivation of the self-energy 

There are several ways of evaluating the self-energy, either by using Wick's theorem (see, e.g. 
Ref. [36]) or by functional derivative. We follow the latter because it is physically appealing 
and give a summary of the steps. For more details we refer to Refs. [47, 59]. We introduce a 
time varying field <j>{x) which functions as a mathematical device and will be set to zero in the 
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end. It is similar to the principle of virtual work in classical mechanics. For a time- varying field 
it is convenient to work in the interaction (Dirac) representation: 

\i> D (r,t)) = U(t,t )\*!>D(r,t )) (40) 

The time development operator U is given by 



f7(Mo) = T exp f t <It4>(t) 



(41) 



4>(t) = J <Pr 0(r, t) $,(r, t)$ D (r, t) (42) 
The relationship between operators in the Heisenberg and interaction representation is 

ip H (x) = UHt,0W D (x)U(t,0) (43) 
The field operator ipo{ x ) satisfies 

i^ D (x) = [j>(x),H(i = 0)] (44) 

so it is the same as the unperturbed (</> = 0) Heisenberg operator. The Green function can now 
be written as 

G( , = . (N°\T[U(^-oo)i D (x)^ D (x')\N ) 
y ' (N°\U(oo,-oo)}\N°) V 

and | A'' ) is the unperturbed groundstate (the interacting groundstate before the application of 
(f>). By taking functional derivative of G with respect to <j> it can be shown that 

SG(i,2) 
6<t>{3) 

Using this result in Eq. (15) we get 



G(l,2)G(3,3 + )-G 2 (l,2,3,3 + ) (46) 



J d3M(l,3)G(3,2) = V H (l)G(h2) + ij d3v(l,3)-^Q (47) 

where 

vil^^vQn-rsMfa-h) (48) 

and Vh is the Hartree potential 

V H (x) = J dx'v(x,x')p(x') (49) 
Using the definition of the self-energy £ = M — Vjj and using the identity 

■.£ (e - lG) _0^- G *£! O (50) 

we have 

£(1,2) = -if d3div(lA)G(h3) 6G 6 ^ 2) (51) 
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This expression is exact. The quantity 8G~ l 18$ is related to the dielectric or response function. 
Multiplying both sides of Eq. (18) by G~ x and keeping in mind that h contains the probing 
field (j> we find 

dE(3,2) 



ft?- 1 (3, 2) 



*(3-4) + ^ S > 



(52) 



We can see very clearly the stucture of the self-energy: Including the first term of 8G~ l /8(j) 
results in the well-known Hartree-Fock approximation. Including further the second term, i.e., 
the change in the Hartree potential, in fact leads us to the GWA because, as indicated, 1+5Vh /6<j> 
is the inverse dielectric function so that 

W{1,2) = Jd3v(l,3)e- l (3,2) (53) 

is the screened Coulomb potential. We therefore have 

E GW (l,2) = iG(l,2)W(l,2) (54) 

The last term in 8G~ l /8<fi is the response of the self-energy to the probing field and this consti- 
tutes the vertex corrections. Thus, apart from the Green function, response functions are very 
important quantities in the calculations of the-self-energy. Note that we have not invoked any 
diagrams or many-body perturbation theory in deriving the GWA or the self-energy in general. 

3.3 The response function 

We have seen that the Schwinger functional derivative technique is a very convenient tool in 
Green's function theory. Here we demonstrate again its usefulness by deriving the equation 
satisfied by the charge response function. It is convenient to use a convention where repeated 
index or variable means a summation or an integral. We first define the polarisation function 
which may be thought of as a response function but with respect to the total potential V — 
<f> + V H : 

P(12) = 4^ 
nL,Z) SV(2) 

■ *C(1,1+) 
1 8V{2) 

= -iG(l,3) ^gj 4) g(4,l+) 

= tG(l,3){<5(3-4)J(3-2) + ^^}G(4,l + ) 
= «G(1,3)A(3,4,2)G(4,1 + ) (55) 
where A is known as the vertex function: 

A(l,2,3)=f(l-2)f(l-3) + ^^ (56) 

We have used the identity in Eq. (50) and Eq. (52) for 8G~ 1 /8(f>. From 

_i 8V_ 

€ ~ 8(j> 

= 1 + V H 



- 451 - 



Sp SV 

= l + vPe~ l (57) 

we find 

c _1 ='[1-«P] _1 or e = \-vP (58) 
We now calculate the response function: 

Sp(l) 



12(1,2) = 



*P(1) ^(3) 



*V(3) <ty(2) 

= P(l ) 3)e" 1 (3,2) (59) 



Thus schematically 



R = P[l - vP]~ l or R = [1 - Pv] _l P (60) 

This equation for P is exact. If we use a noninteracting P° we arrive at the well-known time- 
dependent Hartree or the random phase approximation (RPA) equation. The RPA has a simple 
physical interpretation: The change in density due to a perturbing external field is given by 
Sp = P°6(6V e xt + SVh), i.e., the system response to the total field as if it were not interacting. 
When 1 — Pv = we expect a new mode of excitation which can be understood as follows: 
A density fluctuation Sp giving rise to a potential Sv = vSp such that PvSp = PSv = Sp, i.e., 
the potential generated by the density fluctuation induces the same density fluctuation. This 
new mode of excitation is usually referred to as a plasmon. Thus a plasmon is a self-sustaining 
density fluctuation. 

By using the spectral representation for G in Eq. (24, the polarization function can also be 
expressed in terms of its spectral representation: 

P(r , r ., u) . f & + r „ ,„, 

7_oo w — U)' — to Jo uj — w' + 10 

where S is proportional to the imaginary part of P and defined to be anti-symmetric in u: 

S(r,r',u) = — — Im P(r, r', ai)sgn(w) 

= £Mr)W^(r)p*(r')% - e(N + l,j) - e(N - 

(62) 

For a noninteracting system, hi -+ ip? cc , pj -4 V^ nocc and e(N- 1, i) -+ ef c , e(N+l,j) -+ ^ nocc . 
Similarly, the response function R can be written in the form of Eq. (61). 

In terms of the many-body states, |M), the exact time-ordered response function can be 
expressed as 



fl(l,2; W ) = £ 
M 



(0|/5(1)|M)(M| / 5(2)|0) (0|p(2)|M)(M| / 5(l)|0) 



uj - Em + Eo + i6 u + Em — Eq - i6 
Thus the poles of the response function yield the exact excitation energies. 



(63) 
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4 The G W approximation 

The GWA is usually regarded as the first term in the expansion of the self-energy in the screened 
interaction W. This point ov view is based on many-body perturbation theory. This is, however, 
not very useful from physical point of view. If we do include higher order diagrams, we do not 
necessarily get better results (in the sense of closer agreement with experiment). Straightforward 
higher order diagrams can give a self-energy with wrong analytic properties which in turns results 
in an unphysically negative spectral function (density of states) for some energies. 

Eqs. (51) and (52 on the other hand are exact and the GWA has been obtained by neglecting 
the vertex 8T,/S<j) in SG~ l /S^>. It has not been obtained by summing diagrams as is normally 
done in many-body perturbation theor}'. The GW self-energy has the same form as the exchange 
operator in the HFA with the bare Coulomb potential substituted by a screened Coulomb in- 
teraction. It is more physical and useful to regard the GWA as a Hartree-Fock approximation 
with a frequency-dependent screening which cures the most serious deficiency of the HFA. 

The GWA has been applied with success to many systems ranging from alkali metals [99, 
77, 82, 83], semiconductors [55, 41], transition metals [8, 11], metal surfaces [30, 31] to clusters 
[85]. The list is by no means exhaustive. 

4.1 Explicit expression for Sgw 

Fourier transforming Eq. (54) we obtain 

£(r,r',u/) - ^-/dw'G(r,r',w + w')W(r,r',w') 

27T J 

= I dJ G(r, r', u + J)v{r - r') 

ZTT J 

+£- [ du' G(r, r', u + u')W c (r, r', J) (64) 
2n J 

The first term gives the bare exchange 

E*(r,r') = -2>(r)/*(r>(r-r') (65) 

i 

and the second term gives the correlation part of the self-energy which has the following spectral 
representation: 

J -co uf — ur — id Jfj, u) — u/ + to 

The spectral representation of the correlation part of the screened Coulomb interaction W c — 
W - v is 

r,r , , ^ f° j i D(r,r',u>') f°° , , JD(r,r',w') 

W c (r,r» = / dJ / ' + / dJ v ' , ' .' (67) 

J-OO U) ~ U)' — id Jo u> — U)' + 10 

Using the spectral representations of G and W c in Eq. ( 64) gives (Appendix B) 

/•oo 

T(r, r';u;) = -sgn(u; - fi) / dio'0{fj, — oj - u>')A(r, r', lj + o/)Z>(r, r', u) 

Jo 

rOO 

+sgn(u> - fi) / dw'0(-(j, + u) — uj')A(r, r',(j - w')D(r,r' ,w) 
Jo 

= -sgn(u,-/x) X)/ i (r)/;(r , )D(r ) r' > ai(o;))(9[a i ( W )] 

i 

+sgn(u, - fi) 5>(r)tf(r') D(r, r', AM) 0[&M] (68) 
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where 



ai(w) = -co + n-e{N (69) 
= w-/*-e(JV + l,») (70) 

The real part of the self-energy can be obtained from the Hilbert transform in Eq. ( 66). 
4.2 Y, G w with a no ninter acting Go 

In practice, most GW calculations have been performed using a noninteracting Green's function. 
The issue of self-consistency will be discussed at a later section. For a non-interacting system, 
we have 

e{N - 1, i) = E{N - 1, i) - E{N - 1) 
= -Si + E(N) - E(N - 1) 

= Si + fJ. (71) 
e(N + l,j) = Ej-n (72) 

fi — i>u with £j < (j, (73) 
gj = ip*, with Ej > fi (74) 

where is the single-particle eigenvalue of the the single-particle state tpi of the non-interacting 
system. The spectral function of Go is 

A°(r,r',u,) = X>(r)</>*(rX" - e l ) (75) 

i 

since A is real. For a non-interacting system, the spectral function consists of delta function 
peaks at the single-particle eigenvalues whereas in the interacting system, apart from the quasi- 
particle peak, the spectral function may have a satellite structure. 
The spectral function of P for the non-interacting system is 

5°(r, v>,u) = ]T £ ^(r)^V)^(r)V;(r') 6[u - (ej - *)] (76) 

i</ij>/i 

and the spectral function for Hgw is (Appendix B) 

' Ei< M M^ii^Dir, r', ei - uj)0{Ei - u) for u < \l 

< 

. T,i >fl Tpi{r)Mr')D{r, r', u - ej)0(u> - a) for u > \t 
The bare exchange part becomes 



r(r,r» = - 



(77) 



E*(r,r') = -X^(r)V;(r')^(r-r') (78) 

i</i 

4.3 Coulomb hole and screened exchange (COHSEX) 

A physically appealing way of expressing the self-energy is by dividing it into a screened-exchange 
term Ssex and a Coulomb-hole term Ecoh (COHSEX) [46, 47]. It is straightforward to verify 
that the real part of the self-energy can be written as 

occ 

ReSsEX (r.r'.w) - - £^(r)V*(r')rW (r.r'.w - si) (79) 
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ReScoH (r,r» - Y]^(r)^(r')P f°° daj'^^j (80) 

^ 70 W - £i - Ct/ 

The physical interpretation of Ecoh becomes clear in the static approximation due to Hedin 
(1965a). If we are interested in a state with energy E close to the Fermi level, the matrix element 
(V'lReEcoH(-E')l^') picks up most of its weight from states with energies Si close to E in energy. 
We may then assume that E — e% is small compared to the main excitation energy of D, which 
is at the plasmon energy. If we set E — e{ = 0, we get 

ReScoH (r, r') = \s (r - r') W c (r, r', 0) (81) 

This is simply the interaction energy of the quasiparticle with the induced potential due to the 
screening of the electrons around the quasiparticle. The factor of 1/2 arises from the adiabatic 
growth of the interaction. In this static COHSEX approximation, £coh becomes local. 

4.4 Generalization of the HFA 

The GWA may be regarded as a generalization of the HFA. ^From Eq. ( 64) and using the 
spectral representation of Go, the correlated part of the self-energy can be written as 

£c(r,r';w) - ^^(r)^(r')W c -(ry; W - ei ) 

i<fi 

+ J^WVaOW+try^-Ci) (82) 

i>H 

where 

In short we can write 

S(r,r';a;) = -^V'i(r)V'r(r')M / o(r,r';a;-e i ) (84) 

where 

Wbfcr'jw-ei) = {v(r-r')-W-(ry;u-€i)} 6(f,-ti) 

-W+ir^-u-eJQiti-n) (85) 

Thus, the self-energy in the GWA has the same form as in the HFA except that it depends on 
energy and contains a term which depends on unoccupied states as a consequence of correlation 
effects. Thus the GWA can be interpreted as a generalization of the Hartree-Fock approximation 
(HFA) with a potential Wo which contains dynamical screening of the Coulomb potential. Note, 
however, that Wq is not the same as the dynamically screened potential W . 

As previously mentioned, the HFA does not take into account the effects of screening which, 
for insulators, results in too large band gaps. It can be shown that the GWA gives the correct 
band gap, at least for localized states which are well isolated from the other states. Consider 
the correlated part of the self-energy for an occupied core- like state \p d - 

(^c(e d Md) = (^1^(0)1^) 
occ 

unocc 

+ ^2{Mi\W+(€ d -€i)\i,^ d ) (86) 



% 

27 



oo w + u>' ± iS 



(83) 
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Strictly speaking, the self-energy should be evaluated at the new energy E d = €4+ self-energy 
correction and it is understood to be the case here. If ip^ is localized and well separated in 
energy from other states, then the first term is evidently much larger than the rest. Thus we 
may make the following approximation: 

(1>d\Vc(cdMd)- * (i>jl>d\w-(o)\Md) 

= -\(Md\W e {0)\i>d^d) (87) 

The last step is shown in Appendix A. This is a correction due to the work done on the electron 
by the polarization field from zero to W C (Q) [46, 47]. A similar result, 

+\{Md\wM\Md) 

is obtained for an unoccupied core-like state of the same character so that the energy separation 
of the states is 

A = 4 F - ef F + {Md\^M\Md) 

= (Md\v\i'di>d) + {if>dTl>d\w c {o)\ipdti>d) 

= (MdWrn^d) (88) 

which agrees with the intuitive result that the "gap" is given by the screened Coulomb inter- 
action: A = U « W(0). Since W c (0) is negative, we see that the self-energy correction to the 
HFA raises an occupied state and lowers an unoccupied state. This is still true also for states 
which are not so localized. 

5 Vertex corrections: Beyond GW 

The G WA has proven to be very successful for describing quasiparticle energies for sp and even 
3d systems. Its description for satellite structures, however, is less satisfactory. A number of 
cases which reveal the shortcomings of the G WA in describing the satellite structure are 

• The valence as well as the core photoemission spectra of the alkalis show a series of plasmon 
satellites which are located at multiple of plasmon energies below the main quasiparticle 
peak. The G WA gives only one plasmon satellite and its position is typically 1.5 plasmon 
energy below the quasiparticle peak [48]. 

• The valence photoemission spectrum of Ni shows the presence of a satellite at 6 eV which is 
not obtained within the GWA. Similarly, the GWA appears to be insufficient for describing 
the satellite structure in transition metal oxides and the situation becomes worse in the 
f-systems. 

The presence of only one plasmon satellite in the GWA follows from the fact that the self- 
energy is of first order in W which contains one plasmon excitation through the RPA dielectric 
function. From the diagram, it is clear that a hole or an electron interacts with its surrounding 
by exchanging a plasmon. The too large plasmon energy can be understood qualitatively as an 
average of the first and the second plasmon energies since they carry most of the satellite weight. 
Because the peak in Im S c is at one plasmon energy below the quasiparticle peak, the plasmon 
satellite in the spectral function ends up at «1.5 plasmon energy To get the correct satellite 
energy the peak in Im E c should therefore lie somewhat closer to the quasiparticle energy. 

The problem with the satellite structure in Ni or transition metal oxides is of a different na- 
ture since the satellite energy is much lower than the plasmon energy. In the atomic picture, the 



ground state of Ni is a mixture of the configurations 36^4$ and 3d 8 4s 2 . The final configurations 
after photoemission are 3d 8 4s and 3d 7 s 2 . The former corresponds to the main line (quasiparticle) 
and the latter to the satellite. The two configurations are separated in energy by 6 eV which is 
essentially the Coulomb energy of two d holes. This value is reduced by screening in the solid. 
The presence of satellite at a certain energy implies that the imaginary part of the self-energy 
should exhibit a peak at an energy slightly lower then the satellite energy. Such a two-hole 
excitation state is partly described by the GWA to second order in the bare interaction but the 
GWA mainly describes the coupling of the electrons to the plasmon excitation. A secondary 
plasmon-like excitation is possible within the RPA when the band structure gives rise to two 
well-defined peaks in Im e -1 . 

It is a major challenge to develop a theory beyond the GWA for real systems which can 
overcome the problems described above. Here we describe those approaches which are based 
on systematic diagrammatic expansions. One approach which has been applied with success 
to the alkalis is the cumulant expansion [66, 19, 49]. As described in the next section, this 
approach is suitable for dealing with systems which can be mapped into systems of electrons 
coupled to bosons (e.g. plasmons) where long-range correlations dominate. However, short- 
range correlations arising from multiple-hole interactions on the same site should be better 
treated within the T-matrix approach, described in a later section. In this contribution, we are 
concerned with vertex corrections related to the satellite description rather than quasiparticle 
description. 

The effects of vertex corrections on the quasiparticles have been studied by a number of 
authors. These include a direct calculation of the second-order diagram [27, 28, 21] and vertex 
correction based on the LDA V xc with application to Si [32]. Vertex corrections in the electron 
gas [77, 37] are discussed in a review article by Mahan [79] where it is emphasized that it is 
important to include vertex corrections in both the response function and the self-energy in a 
consistent way [107, 14, 15, 90]. 

5.1 The cumulant expansion 

One of the first applications of the cumulant expansion method was in the problems of X-ray 
spectra of core-electrons in metals [84, 66]. The problem is modelled by a Hamiltonian consisting 
of a core electron interacting with electron-hole excitations and a set of plasmons: 

H = eJc + ]T ekacj^cior + E w q 6 q 6 q 

k<r q 

+ E ^kk'cL^'Cct + £ cct ffq (6 q + b{) (89) 
kk'cr q 

where c is the annihilation operator for the core electron, cj^ is the creation operator for a 
conduction electron of wave vector k, spin a, and energy and 6 q is the creation operator 
for a plasmon of wave vector q and energy w q . The last two terms are the coupling of the core 
electron to the conduction electrons (electron-hole excitations) and to the plasmons respectively. 
The model Hamiltonian without the last term [75, 76] was solved exactly by Nozieres and de 
Dominicis [84] and the full Hamiltonian was solved exactly by Langreth [66] who also showed 
that the cumulant expansion also gives the exact solution. We consider the case with coupling 
only to the plasmon field [72, 73] but not to the conduction electrons. The exact solution is 
given by [66] 

00 e~ a a n 

A± (w) = ]T — 5 (u - e - Ae t nu p ) (90) 

where + refers to absorption spectrum and — to emission spectrum, a = Qq/^p an( l = aw p 
is the shift in core energy due to the interaction with the plasmon field. The spectra consist 
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therefore of the main quasiparticle peak at oj — e + Ae and a series of plasmon excitations at 
multiples of the plasmon energy below the quasiparticle peak. It has been assumed that the 
plasmon excitations have no dispersion although this assumption is not necessary. 

The physics of the cumulant expansion when applied to valence electrons was discussed in 
detail by Hedin [49, 1]. More recently, the cumulant expansion was calculated to higher order 
for a model Hamiltonian by Gunnarsson, Meden, and Schonhammer [43]. 



5.1.1 Theory 

In the cumulant expansion approach, the Green function for the hole (t < 0) is written as 

• G(k,t) = i8(-t)(N\cl(0)c k (t)\N) 

= ie{-t)e- iekt + ch ( k V (91) 

where k denotes all possible quantum labels. For u> < only the first term in equation (26) 
contributes and the hole spectral function is 

A(k,u)<fi) = — ImG (fc, u>) 

7T 

1 r°° . ■ 
= ^J^dte^(N\c{(0)c k (t)\N) 

= -Imi f° dte iut e- iEkt+c ^ k ' t) (92) 
C h (k, t) is called the cumulant. Expanding the exponential in powers of the cumulant we get 



G(M) = G (k,t) 



l + C h (k,t) + ^[C h (k,t)} 2 + 



(93) 



where Go {k, t) — i9(—t) exp(-iekt) . In terms of the self-energy, the Green function for the hole 
can be expanded as 

G = Go + G EGo + GoSGoSGo + . . . (94) 
We may group the cumulant into terms labelled by the order of the interaction n: 



oo 



n=l 



The above equations for G may now be equated and the cumulant can be obtained by equating 
terms of the same order in the interaction. Thus to lowest order in the screened interaction W, 
the cumulant is obtained by equating 

G C h = G EG = AG (1) (96) 

where E — Egvv — iGqW. If Go corresponds to, e.g. the LDA G, then E — Egw — V xc . The 
explicit form of AG' 1 ) is 

AG (1) (1,2) = y"d3d4Go(l,3)E(3,4)G (4,2) (97) 

Using 

occ 

Go (1,2) = i'EtnMKMe-^-VOfa-h) 

n 

unocc 

-i E <t>n{n)K{v2)z- i£n{u - t2] e{h -h) (98) 
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and taking matrix element with respect to an occupied state <pk yields 

r<2 



roo ri2 

AG^(k,t l -t 2 ) = -e- l ^ ti ' t2 n dt 3 dU^^-^ik.h-U) (99) 

Jt\ J -oo 



Without loss of generality we can set £2 = and t\ ■ = t < 0: 

roo roo 
dt 3 
It Jt 3 

The first-order cumulant is therefore 

roo roo 



roo roo 

AG W (M) = -e~ lEkt / dt 3 / dr e !£ * T E(fc,r) (100) 



/oo roo 
dt 1 J dTe t£kT Z(k,T) 

/0 roo 
dt' J dre i£kT X(k,T) + C h {k,0) (101) 

C h (fc,0) is a constant which contributes to an asymmetric line shape of the quasiparticle. 

The cumulant can be expressed as an integral over frequency by defining a Fourier transform 

S(*,t) = J^e-'^Sik,*) 

J Z7T U-00 U} — U}' — tO J^, U) — U)' + ld) 

= iO(-T)e ST r dJe~ iw ' T Y{k,J) 
J —00 

-i0( T )e- 5T duj'e- lu)T T{k,oj') (102) 

The second line uses the spectral representation of £ in equation (66) and the last line is obtained 
from 



du e~ iUT _,,,, T f du" e ~ iu > T 

-iS ' 

= i9(-T)e- iw ' T e ST (103) 



/ — e l . u — u — u 

J 2n u - u> - iS J 2tt u" 



Similarly for the other integral. The self-energy becomes 

/li roo 
dwe _ " T r(A; I w)-t0(T)e-* T / due'^T {k,u) (104) 

The cumulant can now be expressed as, keeping in mind that t < 0, 

cft'y dre ie " r i;(/t,r) 

rf** jf dre iekT \i9{-T)e 6r j due~ iuJT T{k,u) 



-09 (r) e~ 6T J°° du e~ luJT T {k, u) J 

- /°dt' Tdr ^ dwe i(£ *-'"- ia ) T r(Jfe,w) 
/■0 /-oo /-oo 

+ dt' dr due l{6k -" + >V T T(k,u>) 
h Jo Jfi 
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dui — — r (fe, oj) 

oo Ek-OJ -10 



doJ- 

y. Ek-OJ + 10 



u l dt 'l 

•M/: 

ry. e *(e* -w-t<$)t _ i 
+ / ^2" 



■fa r < k - u \, + r*a r <*•"».,) 

oo £jfc - W - ZO A £fc - W + W J 



oo (e^ — oj — iS) 
dY, h {k,uj) 



T(k,oj) 



-iE{k,s k )t + 



f 

J -c 



+ / doj- 



du> 

— ui— iS)t 



oo {ef. — u) — isy 



L)=£ k 

T{k,u>) 



The last line is obtained from 



f y ^ , 

V-oo (g/fc — oj' — iS) 



doj J-qo (oj — oj' — iS) 
dZ h {k,u>) 



U=£ k 



doj 



U=£ k 



where 



iS 



7-00 w - w' - z< 
We can also evaluate C ft (fc, 0) : 

/•oo roo ( rfi 

C h (k,0) = i dt'l dre i£kT {id(-T)e dT doje- iuJT T{k,oj) 

J0 Jt' I J-<x 

-%9 (r) e-* T doj <T iur T (A, w) J 

/•oo /-oo /-oo 

= / dt'l dr doje^ ek - w+lS)T T(k,oj) 

J0 it' Ja 



roc r( 

= i / dt' 
Jo Jn 

- -f 



gL>- — T (A;, w) 

6fc - oj + to 



".fa- r^" 1 



(e fc - a; + z<5)" 



dE? (A, a;) | 
5a; 



where 



°°doj' T(KJ) 



(105) 



(106) 



(107) 



(108) 



(109) 



oj — oj' + iS 

The total cumulant can be conveniently divided into a quasiparticle part which is linear in t and 
a satellite part: 

C h = C% P + C% (110) 

Thus collecting the linear terms in t from equations (105) and (108) we get the cumulant con- 
tribution to the quasiparticle: 



C% P {k, t) = {ia k + 7fc ) + (-iAe fc + m)t 



(111) 
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where 

dS(fc,w) 
«afc + Ik = ' 
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, Ae fc = ReE(k,e k ), % = Im£(k,e k ) (112) 



7§(M)=/ do;- ^(fc.u;) (113) 



The contribution to the satellite is 

e i(e k -u-i5)t 

oo (ejt - w - iS) 
A similar derivation can be carried out for the particle Green function 

G (k, t > 0) = -id (t) e ~WCP(k,t) (n4) 

where k labels an unoccupied state. The result is 

C QP (k,t) = {ia k + lk ) + {-iAe k - rj k )t (115) 

poo J(e k -u+i6)t 
C p s (k,t)= duj- — -jr(fc,w) (116) 

J li (E k — U + 10) 

It is physically appealing to extract the quasiparticle part from the Green function: 

G h QP {k,t) = ^(-Oe-**' 4 ^*' 

= ^(-<)e iQfc+7fc e(-^+^) f , ^^ejt + Aefc (117) 

The spectral function for this quasiparticle can be calculated analytically: 

a tu ^ \ e ~ 7 * r) k cos a k - (uj - E k ) sin a k 

A QP (k,uJ<Li) = — 2 2 (118) 

Thus we can see that the quasiparticle peak is essentially determined by the GW value. 
^From equation (91) we have for t < 

<JV|4(0)cjt(t)|JV) = e -^ +c "W) (119) 

By analytical continuation to t > and using equation (92) the spectral function can be rewritten 
as 

A (fc, w) = -!- r dte iu,t e- iEkt+c ^ k ' t) (120) 

Z7T J oo 

where C h for positive t is obtained from C fc *(— t) = C h {t) since A(/c,o;) must be real. 

The total spectra can be written as a sum of Aqp and a convolution between the quasiparticle 
and the satellite part: 

A (k t u>) = Aqp (k,u) + ~ r dte^-^^ "^ 1 \e c s^ - l] 

2ir J-oo J 

= Aqp (k,u) + Aqp (k,w) * A s (fc,w) (121) 

where 

A s (k,u) = ~ j 'dte^^s^ -l} 

= ^/ctte^[cfg(M) + ^[cg(fc,0] a + ...} (122) 
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The second term Aqp * As is responsible for the satellite structure. The Fourier transform of 
Cg can be done analytically 

/OO fff 
^_j»t c h^ t) (123) 

= J^{e^Cl(k,t)+e-^C h s (k,-t)} 



-oo 

f° fit rt 1 „i(e k +u>-oj'-iS)t 

/ ?- d "'~ *r(k,u')+c.c. 



(124) 



Integrating over t gives 

C h s (k,u<0) = ±Imf dJ T } k,U},) 

7T J-oo (si. - UJ 1 — iS) (Sir + i 



-oo (ek - u' — iS) (sk + u — u' - iS) 



(u/' — cj' — i£) (efc + w — w' — i<$) 

5 r r (A, e fc + u) r (it, a/') 



) ( T(k,e k +u) | r(fc,o/') 1 

Dw" \ U 1 ' — U - £k £ k + U) — U}" ) 

r (A;, e k + u) - r (*, - wr (A;, e*) 



w 2 



(125) 



As follows from equations 118 and 121, the quasiparticle energy in the cumulant expansion 
is essentially determined by E k , which is the quasiparticle energy in the GWA. 
Let us apply the cumulant method to the Hamiltonian 

H = eJc + uj p tfb + gcc* (fet + b) (126) 

which is a simplified version of the model Hamiltonian discussed previously. First we must 
calculate the self-energy to first order in the plasmon propagator 

U — Up + 10 U! + Ulp — 10 

The effective interaction between the core electron and its surrounding is g 2 D (see, e.g., Inkson 
[59] ) . The self-energy of the core electron as a result of the coupling to the plasmon is then 

£ (w) ■ = 4~ I ( w + w ') 9 2 ° W) 
Lit J 

If^ — J_/_J_^_^Jl 

In J u + u — e — id \u' — u p + id u ' + u p — id J 



27 

9 2 



(128) 



u + Up — e — iS 
The spectral function for £ is given by 

r (w) - g 2 5 {u - e + Up) (129) 
and the self-energy correction to the core eigenvalue is 

Q 2 

As = — (130) 

Up 
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The derivative of £ at u> = e is 



C' 1 (Jfe, 0) = and C% (k, t) becomes 







- 









u)— iS)t 

-oo (e — U3 — iS) 
The spectral function is thus given by 

1. . f° 

-00 



2 r (w) = ( ^- I e*""' 



TT J-oo 

n J-oo n n! 



— e 



^ y> 1 / 9 
=o 

which is precisely the exact solution [66] 



00 i / \ 2n 



(131) 



(132) 



(133) 



5.1.2 Comparison between the cumulant expansion and the G WA 

To identify vertex corrections contained in the cumulant expansion, we compare it with the 
G WA. Direct comparison in the self-energy diagrams is, however, difficult if not impossible. 
This is because the cumulant expansion is an expansion in the Green function, rather than 
in the self-energy. It is therefore more appropriate to compare the Green functions in the 
two approximations. In figure 10 the Green function diagrams are shown to second order in 
the screened interaction, which should be sufficient for our purpose. The cumulant expansion 
diagrams are obtained by considering the three possible time-orderings of the integration time 
variables t' in C 2 (k,t) with C(k,t) given by equation (122). As can be seen in the figure, 
the cumulant expansion contains second-order diagrams which are not included in the GWA. It 
is these additional diagrams that give rise to the second plasmon satellite and they are quite 
distinct from the second-order diagram common to both approximations. The interpretation of 
the latter diagram is that a hole emits a plasmon which is reabsorbed at a later time and the 
hole returns to its original state before plasmon emission. This process is repeated once at a 
later time. Thus there is only one plasmon coupled to the hole at one time. In contrast, the 
other two diagrams, not contained in the GWA, describe an additional plasmon emission before 
the first one is reabsorbed, giving two plasmons coupled to the hole simultaneously. Similar 
consideration can be extended to the higher-order diagrams. 

The cumulant expansion contains only boson-type diagrams describing emission and reab- 
sorbtion of plasmons but it does not contain diagrams corresponding to interaction between 
a hole and particle-hole pairs. This type of interaction is described by the ladder diagrams. 
For this reason, the cumulant expansion primarily corrects the satellite description whereas the 
quasiparticle energies are to a large extent determined by the GWA as mentioned before. 



5.1.3 Applications 

The cumulant expansion was applied recently to calculate the photoemission spectra in Na and 
Al [13]. The experimental spectra consist of a quasiparticle peak with a set of plasmon satellites 
separated from the quasiparticle by multiples of the plasmon energy (figure 11a). The spectra 
in the GWA shows only one plasmon satellite located at a too high energy, approximately 1.5 
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Up below the quasiparticle (lib)) which is similar to the core electron case. The cumulant 
expansion method remedies this problem and yields spectra in good agreement with experiment 
regarding the position of the satellites as can be seen in Fig. lib. The relative intensities of the 
satellites with respect to that of the quasiparticle are still in discrepancy. This is likely due to 
extrinsic effects corresponding to the interaction of the photoemitted electron with the bulk and 
the surface on its way out of the solid resulting in energy loss. These are not taken into account 
in the sudden approximation. 

When applied to valence electrons with band dispersions the cumulant expansion does not 
yield the exact result anymore as in the core electron case. Surprisingly, the numerical re- 
sults show that the cumulant expansion works well even in Al with a band width of mil eV. 
Considering its simplicity, it is a promising approach for describing plasmon satellites. 

5.2 T-matrix approximation 

In many strongly correlated systems, such as transition metal oxides, the photoemission spectra 
often show a satellite structure a few eV below the main peak. The origin of this satellite is 
different from that of the plasmon-related satellite which is also found in sp-systems and which 
usually has a much higher energy. The additional satellite found in strongly correlated systems 
is due to the presence of two or more holes in a narrow band after a photoelectron is emitted, 
i.e. it is due to short-range rather than long-range correlations. An illustrative example is 
provided by Cu and Ni. In Cu, the 3cf-band is fully occupied so that after photoemission there 
is only one d-hole in the system corresponding to the 3d 9 configuration resulting in no hole-hole 
correlation. Consequently there is no satellite either and a single-particle theory is sufficient 
to describe the electronic structure of Cu. On the other hand, the ground state of Ni already 
contains a configuration with one d-hole so that after photoemission, the final state contains a 
configuration with two d holes. Since the two holes are localized on the same atomic site, there 
will be a strong d-d interaction resulting in the well-known 6 eV satellite which in the atomic 
picture corresponds to the energy of the 3d 8 configuration. 

So far there is no good ab initio scheme for dealing with short-range correlations. Most works 
have been based on model Hamiltonians which have given important insights into the underlying 
physics but which contain adjustable parameters, preventing direct quantitative comparison with 
experiment. The GWA takes into account long-range correlations through the RPA screening 
which determines to a large extent the quasiparticle energies. It is known that the GWA works 
well for quasiparticle energies but from a number of calculations [8, 11, 12] it is clear that the 
GWA has shortcomings in describing satellite structures. Even the plasmon-related satellites 
are not well described as discussed in the previous section. A natural extension of the GWA 
is to include short-range hole-hole correlations. This type of interactions seems to be suitably 
described by the T-matrix approximation [60]. Previous calculations based on the Hubbard 
model showed qualitatively that the T-matrix theory was capable of yielding a satellite structure 
in Ni [69, 70, 86]. Most other works in Ni have also been based on model Hamiltonians [104, 
98, 24, 58]. T-matrix calculation on a two-dimensional Hubbard model is also found to improve 
the satellite description of the GWA [105]. In this section we develop a T-matrix theory for 
performing ab initio calculations on real systems [96, 97]. The method is applied to calculate 
the spectral function of Ni for which 6 eV satellite is not obtained in the GWA [8]. 

An extension of the T-matrix theory including Faddeev's three-body interaction [33] was 
made by Igarashi [56, 57] and by Calandra and Manghi [24, 80]. The theory has been applied 
to Ni [58] and NiO [80] within the Hubbard model. 
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5.2.1 T-matrix 

The T-matrix is defined by the following Bethe-Salpeter equation: 

T = U+UKT 

= [1 - UK]~ l U (134) 

where U is a screened Coulomb interaction and K is a two-particle propagator. Diagrammati- 
cally, the multiple scattering processes are shown in figure 12. It is evident from the figure that 
the T-matrix describes multiple scattering between two holes or electrons. 
The full expression for equation (134) is 

T ffff /(1,2|3,4) = £7(1 - 2)5(1 - 3)6(2 - 4) 

+ U(1 -2) J dl^2'/W(l,2|^2')T^(l',2'|3,4) (135) 

where we have used a short- hand notation 1 = (ri,*i) and a labels the spin. The kernel K or 
the two-particle propagator is given by 

ff^(l,2|l',2') = iG {\',l)G ff ,{2\2) (136) 

where G a is a time-ordered single-particle Green function 

tG ff (l',l) = (N\T\^{\')fa{\)\\N) (137) 

We have assumed that U is an instantaneous interaction which means that ti = t 2 , £3 = £1, and 
ty = ty- Without loss of generality, we may set t\ — t 2 = and £3 = i 4 — t. Equation (135) 
then becomes 

T <7£r /(ri,r2|r3,r 4 ;<) — 

U(ri - r 2 )<J(ri - r 3 )S(r 2 - r 4 )<$(£) 

+ U{r x -r 2 ) J d 3 r[d 3 r' 2 J dt' K^in^r'ur'^t') 

xT^(r' 1 ,r , 2 |r3,r 4 ;i-0 (138) 

Fourier transformation of the above equation yields 

T (T(T '(ri,r2|r3,r4;u;) = 

U{r\ - r 2 )<5(ri - r 3 )£(r 2 - r 4 ) 

+ CA(n -r 2 ) J d 3 r[d 3 r' 2 J K^u^t' uv' 2 ;oj) 

x T„ <r -(r / 1 ,r' 2 |r3,r 4 ;a;) (139) 

where 

K aa ,(r l ,T 2 \r' l ,v' 2 ;u)=i j ^^(rV n; U - u/)G^(r' 2 , r 2 ; (J) (140) 
The Fourier transforms are defined by 

/oo 
dt e-^f{t) (141) 
-00 

^ g e-F( U ) (142) 



Using the spectral representation of G a 

C , M = r + r (143) 

7_oo U) — CJ' — 10 Jfx 10 - U)' +10 

with a non-interacting A c 

A a {r,r>;u) = £ ViaWVSrWw - £«,) (144) 

i 

an explicit form for the kernel K is given by 

^ff<T'(ri,r2|r'i,r'2;w) = 

~f u-e ia - E ja < - iS 

^(r^)^(rihMr' 2 )^V(r 2 ) 



unocc 



+ T — v * fil2£l^ (145) 

^ W - Bur- £,v + 10 



The first term on the right hand side is due to hole-hole scattering and the second to particle- 
particle scattering. This expression is similar to the RPA polarization propagator but the states 
are either both occupied or unoccupied. 

5.2.2 T-matrix self-energy 

The self-energy can be obtained from the T-matrix which consists of the direct term 

S d (4,2) = -*J dld3G(l,3)T(l,2|3,4) (146) 

and the exchange term 

E x (3,2)=? y<ild4G(l,4)T(l,2|3,4) (147) 

Fourier transformation gives 

S^(r 4 ,r 2 ;w) = -* J d 3 nd 3 r 3 

J % £ <Mn, r 3 ; d - wJiWn, r 2 |r 3 , r 4 ; J) (148) 

rr 1 



E£(r 3 ,r 2 ;u;) = i j d z nd 3 r 4 

I ^ G ^ ruT4;J ~ ^) T cArur 2 \r 3 ,r 4 ;uj') (149) 

We note that for the exchange term there is no summation over the spin since exchange between 
particles of opposite spins is zero. 

The spectral representation of T aa > is given by 



J-oo U> - U)' - 10 J 2 u W - UJ' + 10 



I2tt 

where 



QoA<*>) = --Im T^w) sgn(w - 2/u) (151) 

7T 
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The analytic structure of T is determined by K in equation ( 145) which can be seen by con- 
sidering equation (134). The first order term in the T-matrix is the Hartree-Fock term which 
is independent of frequency and therefore it does not influence the analytic structure. The 
T-matrix without the first order term is then given by 

T = UKU+UKT 

= [1 - UK)~ l UKU (152) 

Thus, the analytic structure is determined by K because it can be shown that [1 - UK]' 1 has 
no poles in the first and third quadrants (the poles of K lie in the second and fourth quadrants). 
Using the spectral representations of G and T the self-energy can be written explicitly as 

Im E^(r 4 ,r 2 ;w > n) = 

/occ 
d 3 rid 3 r 3 J2 V'k'n'<r'(ri)V'k'n'a'( r 3) 
k'n'a' 

x Im T^ri, r 2 |r 3 , r 4 ; w + £k'n'<r') + Ek'nV - 2/x) (153) 



Im S ff (r 4 ,r 2 ;w < /i) = 

/unocc 
d 3 nd 3 r 3 VVnv(ri)Vk'n'<T'( r 3) 



k'nV 

x Im T ff / ff (ri,r2|r3,r 4 ;a; + £k'n'<r') 9(-w-£kW + 2/x) (154) 

A similar expression for the exchange part can be easily derived by interchanging r$ and r 4 in 
the T-matrix. The real part of E is obtained from the spectral representation 

J~oo U) — U) — 10 J n LJ — W + id 

where 

r^(w) — Im E^u;) sgn(w — fi) (156) 

The screened potential U is in general frequency dependent. For a narrow band of width 
A, the time scale for a hole to hop from one site to a neighbouring site is determined by 1/A 
which is large and in frequency space it means that the largest contribution to the T-matrix 
comes from u w 0. Physically, the interaction is essentially instantaneous within the same site 
and therefore it is justified to use a static screened interaction. This static approximation will 
be used in the present work. 

The one- electron spectral functions are obtained from 

A , x i IReEjMI 

l[(J) ^|u;-e i -ImE i (a;)|2 + |ReE < (a;)|2 1 } 

with the shorthand notation i = \cna. We have assumed in the above expression that the 
self-energy is diagonal in the LDA wavefunctions which are used to construct G. 



5.2.3 Double counting and the total self-energy 

To calculate the total self-energy, we add the T-matrix self-energy to the G W self-energy. How- 
ever, this leads to double counting since the second order term in the T-matrix is already 
included in the GW self-energy. A straightforward subtraction of the second order term leads 
however to wrong analytic properties of the self-energy and consequently gives some negative 



spectral weight. To solve this double-counting problem, we proceed as follows. We divide the 
bare Coulomb interaction into the short-range screened potential U and a long-range part vl 
((von Barth 1995): 

v = U + v L (158) 
The correlation part of the GW self-energy is schematically given by 

T-gw = GvRv (159) 

where R is the total RPA response function 

R = (l-Pv)- 1 P (160) 

and P is the RPA polarization function. Using the above division of the Coulomb potential we 
obtain 

E C GW = GURv L + Gv L RU + Gv L Rv h + GURU (161) 

The last term is then subtracted out to avoid double counting. It has the same analytic structure 
as the GW self-energy, but since U is smaller than the bare v, this term is always smaller than 
the GW self-energy for all frequencies. This guarantees that the resulting self-energy has the 
correct analytic properties as in the GWA. Numerically this term turns out to be very small. 
Thus, according to the above scheme the total correlated self-energy becomes: 

Z c gwt = *w + Sr - GURU ( 162) 

where is given by 

S^ = S d + S x (163) 
with S rf and S 1 given by equations ( 148) and ( 149) respectively. 

5.2.4 Numerical procedure 

Since the T-matrix describes scattering between localized holes or particles, it is suitable to work 
with basis functions which are also localized such as the linear muffin-tin orbitals (LMTO). In 
the atomic sphere approximation, the LMTO basis functions are of the form (central cell) 

X £ L (r,k) = ^(r) + £ <^, L ,(r) h R , UtRL (k) (164) 

R'L' 

<j> RL is a solution to the Schrodinger equation in atomic sphere R at an energy in the center of 
the L-band region and 4> RLi is the corresponding energy derivative. h R ,y RL is a constant matrix 
depending on the crystal structure as well as the potential. The Bloch states are expanded in 
the LMTO basis. 

Vw(r)-£xk(r,k)6«i(kr«T) (165) 

RL 

Inserting V'kno- m equation ( 145), we obtain 

occ 

^(r^kv^H- E E E E E 

kn.k'n' R1L1 R1L1 R 3 L 3 R4L4 

X Rl Ll (r'u k)x^ 2 L 2 (ri, k)x^ Ls (r' 2 , k')x^*L 4 (r 2 , k') 

x b Rl h (kna)b R2L2 (kna)b R3L3 (1W)^ (k W) 
w - ew - £k'n'<r' - i& 

+ unoccupied term (166) 
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So far the above expression is exact. The largest contribution to K arises from the onsite terms, 
i.e. within the same unit cell. Furthermore, since U is short-range, the largest contribution to 
the T-matrix also arises from the onsite terms. For practical purpose, we then neglect the k 
dependence of x so that 

K {r i, r 2 \r' i,r' 2 ; u) = 

£ <(r\ K*A<*PW; w) <'(r' 2 )^'*(r 2 ) (167) 

where 

^ b a (kne)b* (kne)b 7 (k'n'a')b* s (k'n'a') 
K acr '(af3\~f8; u) = - > 

kn.k'n' W " £k ™ ~ £kW ~ " 

-(-unoccupied term (168) 

We have used a short-hand notation a = RL 

Using this expression for K in equation (139) and taking matrix element with respect to 
<pRiLi( r i) Vjfeitfte) and Vr'sLsM Vft,L 4 ( r 4) we obtain a matrix equation for T 

+ J] U(a0\ V v) K aa ,( W \\p; w) r^(Ap| 7 «; «) (169) 

which can be easily solved by inversion. 

U(a/3\rju) = j dlrf2y> a (l) V > /9 (2)«(l - 2)^(1)^(2) (170) 

We note that although U depends on the spin through the basis, it has no physical spin depen- 
dence. 

Having obtained T, it is straightforward to calculate the self-energy. From equations ( 153) 
and equations ( 154), the matrix element of the self-energy in a Bloch state ipk n(J is given by 

Im E£(kn;w > fi) = 

occ 

£ b a (kna)bl{kna)b y (k'n'a')b* s (k'n'a') 

k'n'o 7 a/3, 7i5 

x Im T a > a {aP\j8; lo + £ k 'nv) + £k'nV - 2/i) (171) 
Im Ei{kn;u < fx) = 

unocc 

k'n'v' a/3,76 

x Im T^ct^M; w + ek'nv) #( + 2/*) (172) 

The real part of E is obtained from equation ( 155). It is straightforward to derive the corre- 
sponding expressions for the exchange part. 
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5.2.5 Application to Ni with a model screened interaction 



As a first application, we use the T-matrix theory to calculate the self-energy of ferromagnetic 
Ni. In the case of Ni, the T-matrix should describe repeated scattering of two 3d holes on the 
same site, corresponding to the atomic 3d 8 configuration. The mechanism for the origin of the 
6 eV satellite is discussed in the previous section on Ni. 
A model potential of the form 



L/(lr-r'l) = erfc , (a|r - r/|) 



„ rM ( 173 ) 
|r — r | 

is used to calculate the T-matrix. For the T-matrix self-energy, the position of the peaks in 
Im £ is rather insensitive to the screening parameter. However, the intensity varies with a 
and consequently the satellite position, which is determined by Re S, is modified accordingly. 
Thus, if we were allowed to adjust the screening parameter a, the position of the satellite could 
be shifted arbitrarily. It is therefore crucial to determine the screened interaction U from a 
parameter- free scheme. According to constrained LDA calculations [42] Udd ~ 5.5 eV [98]. 
Consequently in our calculations we choose a = 1.2 which is equivalent to Udd ~ 5.5 eV. 

The calculated spectral functions are displayed in figure 13 and compared with the GW 
spectra. The main peak « 3 eV below the Fermi level is the quasiparticle peak. Two satellite 
structures originating from the T-matrix self-energy, absent in the GWA, can be observed below 
the main peak and just above the Fermi level. The position of the satellite below the Fermi level 
is, however, somewhat too low. This is likely due to the difficulties in determining the correct 
screened interaction in the T-matrix and the neglect of particle- hole interaction in the present 
scheme. 

A new interesting feature is the presence of a peak structure just above the Fermi energy 
which arises from particle-particle scattering. At first sight, these scattering processes are ex- 
pected to be insignificant, since the number of unoccupied states is small which leads to a small 
T-matrix for positive energies. However, our results point to the importance of matrix element 
and bandstructure effects, usually neglected in the Hubbard model. As may be seen by an 
examination of equation ( 153) there is a sum over occupied states which amplifies the small 
contribution from the T-matrix. It may be possible to measure the satellite structure above 
the Fermi level in an angular resolved inverse photoemission experiment by choosing certain 
k-vectors, where the quasiparticle peak is well separated from the satellite. 

There is a significant difference between the majority and minority self-energy as reflected 
in the spectral function in figure 13. The probability of creating a hole in the majority channel 
is larger than in the minority channel, since the former is fully occupied. The virtual hole can 
mainly be created in the minority channel due to the presence of 3d unoccupied states just above 
the Fermi level. This implies that the majority 3d self-energy will be larger than the minority 
one which reduces the exchange splitting and the bandwidth, improving the GW result. It is 
found that the T-matrix self-energy reduces the exchange splitting by « 0.3 eV [96, 97] for states 
at the bottom of the 3d band and thus also improving the GW bandwidth. But contrary to 
model calculations [69, 70], the T-matrix has a significantly smaller effect on the bandwidth. 



6 Self-consistency and total energies 
6.1 Self-consistency 

The current practice of performing GW calculations is to choose a non-interacting Hamiltonian 
(usually the LDA) and to construct the Green function Go from the eigenstates of the Hamil- 
tonian and to calculate the response function. These are then used to calculate the self-energy 
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which is taken as the final result. This procedure is non-self-consistent since the interacting 
Green function G obtained from the Dyson equation G - Go + GoEoG is not necessarily the 
same as Go- To achieve self-consistency, the interacting Green function should be used to form 
a new polarization function P — -iGG, a new screened interaction W, and a new self-energy 
S. This process is repeated until G obtained from the Dyson equation is the same as G used to 
calculate the self-energy. Self-consistency guarantees that the final result is independent of the 
starting Green function Go. Moreover, a self-consistent GW scheme ensures conservation of par- 
ticle number and energy. This is a consequence of a general theorem due to Baym and Kadanoff 
[14] and Baym [15] who proved that approximations for £ are conserving if E is obtained as a 
functional derivative of an energy functional $ with respect to G : 

Conservation of particle number means that the continuity equation 

-$n(r,t) = V-j(r,t) (175) 
is satisfied when n and j is obtained from the self-consistent Green function. Furthermore 

N = -tr r du!mG{uj) (176) 

TT J-oo 

gives the correct total number of particles. Conservation of energy means that the energy change 
when an external potential is applied to the system is equal to the work done by the system 
against the external potential when calculated using the self-consistent G. 

The first self-consistent calculation GW calculation was probably by [29] for a model quasi- 
one dimensional semiconducting wire. The relevance of this model to real solids is, however, 
unclear. Self-consistent calculations for the electron gas were performed recently by von Barth 
and Holm [106] and by Shirley [91]. The calculations were done for two cases: in the first case the 
screened interaction W is fixed at the RPA level, W = Wq, (calculated using the non- interacting 
Go) and only the Green function is allowed to vary to self-consistency and in the other case both 
G and W are allowed to vary (full self-consistency case) [54]. The results of these studies are 

• The band width is increased from its non-self-consistent value, worsening the agreement 
with experiment (figure 14). 

• The weight of the quasiparticles is increased, reducing the weight in the plasmon satellite. 

• The quasiparticles are slightly narrowed, increasing their life-time. 



The plasmon satellite is broadened and shifted towards the Fermi level (figure 15). In the 
full self-consistent case, the plasmon satellite almost disappears (figure 16). 



The main effects of self-consistency are mainly due to allowing the quasiparticle weight Z to 
vary. The increase in band width is disturbing and can be understood as follows. We consider 
the first case with fixed W = Wq for simplicity. First we note that the GW result for the band 
width after one iteration is close to the free electron one. This means that there is almost a 
complete cancellation between exchange and correlation. After one iteration the quasiparticle 
weight is reduced to typically 0.7 and the rest of the weight goes to the plasmon satellite. The 
new Im S, calculated from G obtained from the first iteration has its weight reduced at low 
energy and increased at high energy compared to the non-self-consistent result (figure 17). This 
is due to the sum-rule [106] 

/OO /-OO 
du\ImZ c (k,w)| = Y, / du\ lmW (q, w ) I ( 177 ) 
-oo „ JO 
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which shows that the left hand side is a constant depending only on the prechosen Wo but 
independent of k and self-consistency. For a state at the Fermi level, self-consistency has little 
effects since Im £ has almost equal weights for the hole (w < /i) and the particle part (u > //) 
which cancel each other when calculating Re £ c , as can be seen in equation (66) and as illustrated 
in figure 18. But for the state at the bottom of the valence band, ImS has most of its weight 
in the hole part (see e.g. figure 6) so that the shifting of the weight in ImE to higher energy 
causes Re £ c to be less positive than its non-self-consistent value. A similar effect is found for 
the exchange part which becomes less negative but because the bare Coulomb interaction has no 
frequency dependence, the renormalization factor has a smaller effect on S* so that the reduction 
in S 1 is less than the reduction in S c . The net effect is then an increase in the band width. 
The shifting of the weight in ImS to higher energy has immediate consequences of increasing 
the life-time and the renormalization weight Z (through a decrease in dReH c /du>\) of of the 
quasiparticles and of broadening the plasmon satellite, compared to the results of one iteration. 

When W is allowed to vary (full self-consistencj') the band width becomes even more widened 
and the plasmon satellite becomes broad and featureless, in contradiction to experiment. The 
quasiparticle weight is increased further [54]. These can be traced back to the disappearance of 
a well-defined plasmon excitation in W. The reason for this is that the quantity P = -iGG no 
longer has a physical meaning of a response function, rather it is an auxiliary quantity needed 
to construct W. Indeed, it does not satisfy the usual /- sum rule. The equation e(q, u) p ) = 
determining the plasmon energy is not satisfied any more since the Green function now always 
has weight around w = uj p . This has the effect of transferring the weight in Im £ even further to 
higher energy with the consequences discussed in the previous paragraph. 

A self-consistent calculation for the electron gas within the cumulant expansion has been 
performed by [53]. The result for the quasiparticle energy is similar to the self-consistent GWA 
since the quasiparticle energy in the cumulant expansion is essentially determined by that of 
the GWA as discussed in the cumulant expansion section. The satellite part is little affected by 
self-consistency. 

6.2 Total energies 

So far, the GWA has been applied mainly to calculate single-particle excitation spectra, but 
it is also possible to calculate the total energy [54, 35]. Holm and von Barth [54] calculated 
the total energy using the self-consistent Green function and self-energy in the Galitskii-Migdal 
formulation [38]: 



where ho is the kinetic energy operator plus a local external potential. The total energy turns out 
to be very accurate in comparison with the quantum Monte Carlo (QMC) results of Ceperley 
and Alder [25]. For r s = 2 and 4 QMC gives 0.004 and -0.155 Rydberg respectively while 
self-consistent GW gives 0.005 and -0.156 Rydberg. This is rather surprising since the GWA 
represents only the first term in the self-energy expansion. The reason for the accurate results 
is not fully understood. Applications to other more realistic systems are necessary to show 
if the results are of general nature. It is probably related to the fact that the self-consistent 
GW scheme is energy conserving and it is partly explained by consideration of the so called 
Luttinger-Ward energy functional [74]: 



which is variational with respect to G and it is stationary when G is equal to the self-consistent 
G which obeys the Dyson equation: 




(178) 




(179) 



— = when G — G + G £G 



(180) 
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= 1? — is the grand canonical potential whose stationary value corresponds to the energy 
(minus /J.N) obtained from the Galitskii-Migdal formula and 3> is an energy functional consisting 
of irreducible energy diagrams. It is not clear, however, why the first order energy diagram 
(giving the GW self-energy upon taking functional derivative with respect to G) appears to 
represent a very good energy functional. Furthermore, the chemical potential calculated from 
\i — dE/dN is in agreement with the value obtained from \i — kp/2 + E(kp,0). The particle 
density n = 2 £ k dujA(k, w) yields n = hp/ (Sir 2 ) , i.e. particle number is conserved, as 
proven by Baym. It can also be shown that with a fixed W = Wo particle number is also 
conserved [54]. The conclusion is that fully self-consistent GW calculations for quasiparticle 
energies should be avoided. It is more fruitful to construct vertex corrections (beyond GW) or 
to perform partially self-consistent GW calculations where the choice of G and W is physically 
motivated. For instance, one can fix W at the RPA level and modify G by using quasiparticle 
energies but keeping the renormalization factor equal to one. Or one could choose a single- 
particle Hamiltonian such that the GW quasiparticle energies are consistent with the single- 
particle eigenvalues. 

Recently, Almbladh and Hindgren [51] discovered that the GW total energies of the electron 
gas calculated using the Luttinger-Ward formalism [74] are in very good agreement with the 
QMC results of Ceperley and Alder and the self-consistent results of Holm and von Barth. Due 
to the variational nature of the Luttinger-Ward functional, the Green function need not be the 
self-consistent one. A generalization of the Luttinger-Ward formalism by treating both G and W 
as independent parameters has also been developed [2]. Using the non- interacting G — Go and 
the plasmon-pole approximation for W gives almost as good results [51] as the self-consistent 
results of Holm and von Barth [54] or the QMC results of Ceperley and Alder [25]. These results 
are very encouraging and applications to real inhomogeneous systems are now in progress to test 
the validity of the scheme. 

7 Computational method 

The calculation of the self-energy involves the calculations of the following quantities: 

1. A self-consistent bandstructure which is the input to the self-energy calculation. In prin- 
ciple we may use any non- interacting Hamiltonian but we use the LDA in practice. 

2. The bare Coulomb matrix v. 

3. The polarization function P 

4. The response function R and the screened Coulomb matrix W. 

5. The self-energy E. 

7.1 Basis functions 

To construct a minimal basis, let us consider the polarization P. Due to the symmetry of P 
with respect to a lattice translation T, 



P(r + T, r' + T, w) = P(r, r', w), 



(181) 



it follows that P may be expanded as follows: 




(182) 



kij 



where Bfc are any basis functions satisfying Bloch's theorem and they are normalized to unity 
in the unit cell with volume ft. All other quantities depending on two space variables such as 
G, W, and S can be expanded in a similar fashion. The polarization matrix Pij is given by 

fly(k, W ) = jA|dV5yr)P(r,r' )W )B k .(r') 

= E / rf3r / dV s L( r + T ) p ( r + T ' r ' + T '' w )%/( r ' + T ') 

= Y. I d * r I e~ ik (T " T,) ^.(r)P(r, r' + T' - T, wJB^V) 

= N 2 [ <Pr [ <Pr' Bi i ir)P(ry,u,)B kj (T>) (183) 
Thus the calculation of P for the crystal is reduced to a unit cell. ^,From Eq. (76) we then get 

It n</in'>/i 

x^-(^k +q) n'-^kJ] (I8 4 ) 

In the above expression, the wave function ^ is normalized to unity in the unit cell. The 
matrix elements are given by 

<V> k+ qn< U-kn^tf) = ^ ^ +q „^ k n^' ( 185 ) 

The problem is to choose a minimal number of basis functions needed to describe the response 
function. It is clear from the above expression that that the basis functions must span the space 
formed by products of the wave functions. 

If the wave functions are expanded in plane waves, as in pseudopotential theory, then the 
basis functions will be products of plane waves which are also plane waves, in which case 

exp[i(k + G) • r] . cO 

VH ' GG' 

One simply needs to have a large enough G vector in order to have a complete basis. The ad- 
vantages of plane-wave basis are that matrix elements can be easily calculated and the Coulomb 
matrix is simple since it is diagonal, ^QQ'(k) = 4w6qq> /\k'+ G| 2 . Other advantages are good 
control over convergence and programming ease. There are, however, serious drawbacks: 

1. It is not feasible to do all electron calculations. In many cases, it is essential to include 
core electrons. For example, the exchange of a 3d valence state with the 3s - 3p core states 
in the late 3d metals is overestimated by the LDA by as much as 1 eV which would lead 
to an error of the same order in the pseudopotential method. 

2. The size of the response matrix becomes prohibitively large for narrow band systems due 
to a large number of plane waves. 

Moreover, the plane waves have no direct physical interpretation. 

To overcome these drawbacks, we use the LMTO method which allows us to treat any system. 
The LMTO method uses a minimal number of basis functions and we carry over the concept 
of minimal basis in bandstructures to the dielectric matrix e. Instead of a planewave basis, we 
use a "product basis" which consists of products of LMTO's. As will be clear later, the product 
basis constitutes a minimal basis for e within the LMTO formalism. A method for inverting 
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the dielectric matrix using localized Wannier orbitals instead of planewaves has also been used 
in the context of local field and excitonic effects in the optical spectrum of covalent crystals 
[44, 45]. 

In the LMTO method within the atomic sphere approximation (ASA), the basis functions 
are given by 

XRLu{r) = <j)RLu{r) + Yl <t>R>L'v'{r)hR<L><s>,RLu (186) 
R'Uu< 

The orbital ^ is a solution to the Schrodinger equation for a given energy e„ and <f) is the energy 
derivative at e v . The label R denotes atom type, L — Im denotes angular momentum, and v 
denotes the principal quantum number when there are more than one orbital per L channel [9]. 
Therefore a product of wave functions consists of products of the orbitals <j> and <f> [10]: 

^RLu^RL'u'i §RLv§RL'v' \ 4>RLv4>RVv' (187) 

This means that these product orbitals form a complete set of basis functions for the polarization 
function and also for the response function and the self-energy as discussed below. ^From Eq. 
(60) we have 

R = P + PvP + PvPvP + ... (188) 
Writing P — \i)Pij{j\-> the second term can be written 

PvP=J2\i}P ij (j\v\k)P kl (l\ (189) 
ij,kl 

Similar expressions can be written down for the other terms and we can therefore write R = 
Y^ij i-e. P and R span the same space. 

The self-energy of a given state Vk n in the GWA schematically has the form 

£(kn,w) = < ^ k JGWh/> kn > 

= < iftto^MV'V'kn > + < V'knV'l^l^kn > (190) 

The two iptp come from the G. We see that E is sandwiched by products of two wave functions 
and it is therefore sufficient to have v expanded in the product orbitals. Thus, the product 
orbitals in Eq. ( 187) form a complete basis for a (7 calculation. 
We reduce the number of product orbitals in three steps: 

1. We neglect terms containing <j> since they are small (0 is typically 10% of <j>). This reduces 
the number of product functions by a factor of 3. In some cases, it may be necessary to 
include </>. 

2. If we are only interested in valence states, then there are no products between conduction 
states. Therefore, in sp systems, products of 4>d<l>d can be neglected and similarly in d 
systems, products of 4>/(f>f may be neglected. 

3. The remaining product orbitals turn out to have a large number of linear dependencies, 
typically 30 — 50%. These linear dependencies can be eliminated systematically, which is 
described below. 

A product orbital is defined by 

bi{r) = 0R Ll »</>Rx V (r) 

= < PTL v ( r )V f Ri>A r )yi>(*)yv(*) ( 191 ) 



- 475 - 



where i — (R, Lis, L'v'). Due to the ASA, this function is only non-zero inside a sphere centered 
on atom R. Thus, there are no products between orbitals centered on different spheres. For a 
periodic system we need a Bloch basis and perform a Bloch sum 

b ki(*) = E e£k T< ^RL,( r - R - T )4>RLv( r - R - T) (192) 
T 

The k dependence is in some sense artificial because the function has no overlap with neigh- 
bouring spheres, similar to core states. After leaving out unnecessary products (step 2), we 
optimize the basis by eliminating linear dependencies (step 3). This is done by orthogonalizing 
the overlap matrix 

Oij =< bi\bj > (193) 

Oz = ez (194) 

and neglecting eigenvectors z with eigenvalues e < tolerance. The resulting orthonormal basis 
is a linear combination of the product functions: 

Bi = Y,bjZj U (195) 

3 

and typically we have ~ 70 - 100 functions per atom with spdf orbitals. The above procedure 
ensures that we have the smallest number of basis functions. Further approximations may be 
introduced to reduce the basis. 

In Table 1 we show a completeness test for the basis. The slight discrepancy for high lying 
states is due to the neglect of 4> in the product basis which become more important for the broad 
high lying conduction states, and also because the optimization procedure puts less weight on 
those products which have smaller overlap. This is not crucial for two reasons: the matrix 
elements become smaller for the higher states, and in relation to GW calculations, there is a 
factor of 1/cj which makes the higher states less important. 

As applications of the product basis, we have calculated as examples the energy loss spectra 
of Ni (Fig. 2), Si (Fig. 3). and NiO (Fig. 8). 



7.2 Special directions 

When calculating matrix elements using the product basis, we encounter angular integrals of 
the form 

J dSl yi^VLiVLiVLA ( 196 ) 

Analytic evaluation of these integrals are computationally expensive. Instead, they are calculated 
by using special directions which are analogous to Gaussian quadratures in one dimension. In 
general, Gaussian integration over a unit sphere means that we try to find M directions f2; and 
weights Wi such that 

M r 

£ WVlW) = -j= for 0<l< Imax, -l<m<l (197) 
<=i v4tt 

Gaussian accuracy is achieved when the number of correctly integrated spherical harmonics is 
equal to the number of free parameters which is 3M — 2, since the spherical harmonics transform 
among themselves under rotation so that one direction can thus be taken to be the z direction 
and the sum of the weights is one. In one dimension, Gaussian accuracy is always achieved 
and the mesh points are uniquely determined. In two dimensions, it is rarely achieved although 
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one usually comes rather close. There are often several sets of directions that yield the same 
accuracy. We have found by minimization of 

M c 

£l£«m(fli)-4=l 2 ( 198 ) 

a set of 62 cubic directions which correctly integrates all spherical harmonics up to and including 
I — 11 plus a few more for the I = 12 harmonics, a total of 168 functions. Full Gaussian accuracy 
would mean 3 x 62 - 2 = 184 functions and the cubic constraint only leads to a 9% less effective 
integration formula. The directions and the corresponding weights are 

Direction W eight / 4ir 



(199) 



(1, 0, 0) 0.130 612 244 897 931 /6 

(1, 1, 0.128 571 428 571 554/8 

(0.846 433 804 070 399, 0.740 816 326 530 515 /48 
0.497 257 813 599 068, 0.190 484 860 662 438) 

plus all possible cubic variations of these (sign changes and permutations). We have also found 
a larger set with 114 cubic directions which integrate up to I — 15. The directions and weights 
are 

Direction W eight [Atx 

(1, 0, 0) 0.076 190 476 192 774 /6 

(1, 1, 0)/y/2 0.137 357 478 197 258 /12 

(0.733 519 276 107 007, 0.344 086 737 167 612 /48 
0.570 839 829 704 020, 0.368 905 625 333 822) 

(0.909 395 474 471 327, 0.442 365 308 442 356 /48 
0.385 850 474 128 732, 0.155 303 839 700 451) 

Using these directions, any angular function can be integrated easily 

dfi/(ft) = £u; i /(fl i ) (201) 



(200) 



7.3 Evaluation of the Coulomb Matrix 

We consider one atom per unit cell for simplicity. Extension to several atoms is straightforward. 
The Coulomb matrix is given by 

where Bqi is normalized to unity in the unit cell. The integrations over the whole space may 
be reduced to integrations over a unit cell Q by using the property 

B qi (r + T) = e i *» ,T Bq i (r) (203) 
and noting that the integration over r' is independent of the origin of r. This gives 

«ii(q) = f d z r f SsB^E^s^B^v) (204) 
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where 



4, e -(q+G)V4a» t(q+G) . (s _ r) 
fi£ (q + G)2 

+ aV e ^ T er/c i Q|s - r : T|) (205) 
^ a|s-r-T| 

The Ewald method has been used to obtain the above decomposition into summations in the 
reciprocal and real space, erfc is the complementary error function equal to (1 — erf), and 
a is an arbitrary constant whose value is chosen to give a fast convergence in the number of 
reciprocal lattice vectors and the number of neighbours. The essence of the Ewald method is to 
add and substract a Gaussian charge distribution which breaks the Coulomb potential from a 
point charge into a short and long range part. The short range part is done in real space and 
the long range part is done in reciprocal space. The main task is to calculate the potential 

•Ms) = / d 3 r Eqis^Bqjir) 

= £pqi(s)*y (206) 

i 

with z given by Eq. (194) and 

Pqi(s) = f^d 3 r bqi(r)E q (s,r) 

= f> qi (s,G) + £> qi ( S ,T) (207) 
G T 



where 



and 



v a (s G) = 4vr e-(q+ G > 2 / 4 " 2 <(q+G ). s f Sr e -i(q+G).r, , , 
PW S ^> ~ (I (q + G) 2 L Cll[ ' 



(208) 



T) = ae* T / d\ ^Ns-r-TI) 
Jn a\s - r - T| 

It is straightforward to calculate Pqj(s, G), since it is a Fourier transform of bqj{r). To calculate 
Pqi(s, T), we use the following expansion formulas 

]^7r^^T-^ ni " )vUi) (210) 

and 

er/(a|s-r|) ^ 4ir , . /„, ,\ 

«V- ri = E 2iTr«(»-.^WwW ( 211 ) 

The coefficients gi{r,s) are determined by numerical integrations using special directions in Eq. 
( 199). 

At the central sphere, bqi(r) has no q dependence and is given by Eq. ( 191), so that 

Pqt(s,T) = ae lC ^Yl w ii{ s T)yL{^) J d9,y L ^y L2 yL (212) 
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where s>p = s — T and 



Wli{ST) = 21 + 1 Jo ^ {p+r _ ^ r ' s T)|Mr) (213) 
We note that the sum over / is cut off by l\ and /2 hi the product function b. Finally, 

«y(q) - / d 3 Sj B qi (s)$ qi (s) (214) 
which is easily done with special directions in Eq. ( 199). 
7.4 Evaluation of the polarization P 

Calculations of P are the most time consuming due to a large number of matrix elements, a 
summation over the Brillouin zone, which is not restricted to the irreducible zone, and a sum 
over occupied and unoccupied states as may be seen from the following expression: 

PifaM = EE E< B q^knlV'k + q,„'>^k + q,n'IV'kn 5 qj) 
k n<nn'>ii 

f 1 1 1 



I W - £ k+q,n' + £ kn + iS " + e k+q,n' " £ kn ~ iS j 

(215) 

For real frequencies, we calculate S° in Eq. ( 76). The 5 function is replaced by a Gaussian 

Six) => eXP " ( ^ /<T)2 (216) 

The self-energy is not sensitive to the choice of a. 

The matrix elements reduce into integrals of four orbitals 

J d 3 r (t>RL i vARL 2 u2 < pRL 3 u 3 4>RL A u i = j dr ^ipju^ipm^ipm^ipRi^ 

x J dn y Ll y L2 y L3 y Li (217) 

The angular integral is calculated by using special directions in Eq. ( 199). 

To obtain the real part of P we calculate the Hilbert transform in Eq. ( 61) using the 
anti-symmetry of 5: 

ReP{u) = f°° du' S<V) 1—^—7 ^—7} (218) 

JO l CO — U>' U) + UJ' ) 

Although the integrand diverges when u' = u, the integral is well-defined because it is a principal 
value integral. In practice, S° is expanded in Taylor series around uj within an interval u — h 
and u> + h. 

For imaginary frequency, we calculate P directly from Eq. ( 215) by setting u> -> iu: 

Pijfaiu) = EE E^q^knlV'k+q.n'XV'k+q.n'IV'kn^qj} 
k n<fin'>n 

~ 2 ( g k+q,n' ~ e kJ 



w2 + ( £ k+q ; n' - £ kn) 2 

(219) 



- 479 - 



Thus P(r, r', w) is real along the imaginary axis although the matrix representation Pij may be 
complex due to the matrix elements. 

The Brillouin zone integration is performed using a simple sampling method. It is also 
possible to use a more accurate tetrahedron method but the replacement of the 6 function by a 
Gaussian is not possible anymore, resulting in a significantly more complicated programming. 

Once we have obtain P, the response function R and the screened Coulomb interaction W 
can be calculated straightforwardly. 

7.5 Evaluation of the self-energy S 

Taking the matrix element of the bare exchange in a Bloch state ^qn we get from Eq. ( 78) 



(220) 



k y 



obtained by expanding the Coulomb potential like in Eq. ( 182). ^From Eq. ( 77) the correlated 
part of Im £ is given by 

' Ek £n'</« Eiy^qnV'k-q.n' l B ki> *W k > w ~ e k-q,n') 

x < 5 lyl^k-q,n'^qn> ~ e k-q,n') /° r u ^ A* 

EkEn'>^Ei>(V'qnV'k-q,„'l- B ki) £ k-q,n' ~ w ) 

x ( B kj#k-q,n'^q") #( e k-q,n' _ w ) for w>n 



TqnM 



(221) 



The real part of S c is obtained from the Hilbert transform (principal value integral) in Eq. 
( 66). As in the case of the polarization, care must be taken when oJ — ui by expanding T in 
Taylor series around u. 

The quasiparticle energy can now be calculated as follows: 



Eqn — £qn + A£ qn (.Eq n ) 

. A r-, / \ . / iTt ,9AEq n (eq n ) 
= £qn + ^^qn(,£qnj + l-C'qn — £<\n) -q^ 



where 



(222) 

AS qn (u;) = (Vqn|Pe S(w) - V* c |Vqn> (223) 

The self-energy correction AS is obtained from first order perturbation theory from Eq. ( 39) 
and the Kohn-Sham equation: 



(H + £)¥ = EV 



(224) 



(Ho + V xc )i> = ej> 



(225) 



where Ho is the kinetic energy plus the Hartree potential. The self-energy correction to £q n is 
given by 



Aeq n — Eqn — £qn 

— Zq n A£q n (e qn ) 



(226) 
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where 



t aASq B (eq w ) 
9a; 



< 1 



(227) 



is the quasiparticle weight. 

The frequency integration of the self-energy may also be performed along the imaginary axis 
[41] plus contributions from the poles of the Green function (Fig. 4). ^From Eq. ( 82) we have 

£qr» = EE E^qnV'k-qyl^ki) ^k^k-q.n'^V) 
k n'<ti ij 



i_ r c 



dJ 



W^(k,a/) 



w + <J - e k _ qn , - iS 



+ EE E^qnV>k-q,n'i 5 ki) ( S kjl^k-q,n^qn) 
k *'>n ij 



1 f°° A ' 

x— / dJ 
2?r L 



oo W + w' - £k-q,n' + 



(228) 



(229) 



We consider the integration along the imaginary axis with a/ — > iu/', a>" real, and along the path 
C (Fig. 1): 



i r°° 

27T 7-oo 



W§(k,o/) 



50 do;" 



/•oo J,." 



+ 



u + iw" - e k _ q n , u - tu" - e k _ q „, 
W£(k,u/) 



+ a/ - £ k-q 1 n' ±i<5 



/•OO 1 

= - / dJ'W^w")- 

JO 7T 1 



W — 



£ k-q,n' 



(^-^k-qy) 2 +-'" 2 
± Wg[k,±(a; - £ k _ qn ,)] 0[±(w - e k _ q>n ,)] *[±(w - a*)] 
* *[±(£k-q,n' - /*)] 



(230) 



The first term is the contribution along the imaginary axis and the second from the poles of G. 
The integrand in the first term is very peaked around u/ = when u> — e k -q n ' 1S sman - To 
handel this problem, we add and substract the following term 

w ~ £ k-g,n> 



= ^A )*— 7 re°' , "" k -1-' ) 'er/c[ Q (o. - £ k ,)] (231) 

z v w fc k-q,n'/ 



When this term is substracted from the integrand in the first term, the resulting integrand is 
smooth and a Gaussian quadrature may be used. 

The GWA has been applied to many systems. Here we show some results for Ni (Fig. 5) 
and NiO (Fig. 7). The self-energy of Ni at the T-point is shown in Fig. 6. The imporatnce of 
the starting Hamiltonian in GW calculations is illustrated in Fig. 9. 
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APPENDIX A 



^From the spectral representation of W c in Eq. (67) we have 

J-oo —{J - id Jo -U' + 10 

f°° , D(ut') 

= -2 / duj'-^-L (232) 
Jo U> — 10 

using the fact that D{ui) is odd. 

H7<oj . f r*« 

2lX J-oo w — 10 



2n J-oo ijj — i 

- j£ 



oo ur — i6 
° , „ £>( W ") , /•«» J „ D(o;'0_- 



^c(O) (233) 



APPENDIX B 



E c (£) = — / G(£? + u)W c {u) (234) 

27T 7_cxd 

The space variables have been dropped out for clarity. Using the spectral representations of G 
and W c in Eqs. ( 24) and ( 67) we get 

A{u>!) 



= rJl d " {/>' 

roo 

+ / dw\ 

J it E + w — w\ -f id 



E + u> — u)\ — iS 



x r du 2 D(u> 2 ) ( : — ^ -1 (235) 

7o I w - W2 + it? u + u> 2 - it] J 

Performing the contour integration in ui yields 

+ r^r^ /^^i, (236) 

The spectral function of S c in Eq. ( 66) is then 

/H roo 
dui du> 2 A(u)i)D(uj2)S{E + U2 -wi) 
-oo Jo 

/•OO /-DO 

+sgn(E - n) du>i du> 2 A{u>i)D(u> 2 )6(-E + u> 2 + u>i) 
Jfi Jo 
roo 

-8gp(E-fi) / duj 2 9{n - E - u 2 )A(E + u> 2 )D{u) 2 ) 
Jo 

roo 

+sgn(E -n) / du 2 6(-n + E- u> 2 )A(E - uj 2 )D{uj 2 ) (237) 
Jo 
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Putting back the space variables and using A in Eq. ( 75) we obtain 

r(r,r',£7) = -sgn(£ - /i) £ *(M ~ £i)^(r)^(r')D(r, r'; £i - £)0te - £) 
+sgn(£ - /x) £ 0( £ , - M )^(r)^(r')D(r, r'; £ - £( -)0(£ - e <) 

i 

(238) 

Finally, 

f Ei< P ^(r)^(r')^(r,r', ei - E)9{e i - E) for E < n 
T(r,r>,E) = \ (239) 
[ Ei>^rW^(r')^(r,r',E- ej )^- £i ) for E > n 
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3s .158114 .158113 .000001 

3p .066833 .066832 .000001 

3p .209174 .209172 .000002 
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.001738 


.000145 


42.29 


.015563 


.015400 


.000163 


42.29 


.018353 
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.006292 


.006040 


.000252 
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.017300 


.016660 


.000640 


73.89 


.011847 


.011409 


.000438 


73.89 


.013877 


.012911 


.000966 



Table 1: The completeness test of the optimal product basis for Nickel. A product of two 
wavefunctions is expanded in the basis: V'^ n ' ! / , ] c ' n / — Ei-^iCi with k — (0 0), ej^ = -1.22 eV 
(the higest valence state at the T point) and k' = (.5 .5 .5). The basis is complete if | £) f |cj| 2 = 
| < V'kJ'^'k'n' > I 2 ' ^ ne numDer of optimal product basis functions is 101 and 82 with and 
without 3s, 3p core states respectively. 
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FIGURE CAPTIONS 



Figure 1: Bandgaps of a semiconductors and insulators calculated within the LDA and the 
GWA, compared with experiment. The data are provided by Dr. Eric Shirley [50]. 



Figure 2: The loss spectra of Ni with (solid line) and without (dotted line) local field compared 
with the experimental spectrum (full circles). Both theoretical spectra are calculated with 4s, 4p, 
3d, 4f, and 5g LMTO orbitals, including an empty sphere at (0.5 0.5 0.5)a and core excitations 



Figure 3: The loss spectra of Si with (solid line) and without (dotted line) local field compared 
with the experimental spectrum (full circles). Both theoretical spectra are calculated with 3s, 
3p, 3d, and 4f LMTO orbitals including core excitations [10]. 



Figure 4: The analytic structure of E c = iGW c for u > Ep (a) and u < Ep (b). Frequency 
integration of the self-energy along the real axis from — oo to oo is equivalent to the integration 
along the imaginary axis including the path C. 



Figure 5: The bandstructure of Ni along FX and TL averaged over the majority and minority 
channels. The solid curves are the experiment and the dotted curves are the LDA (Martensson 
and Nilsson 1984). The filled circles are the quasiparticle energies in the GWA [8]. 



Figure 6: Ni self-energy in the GWA. (a) The real and imaginary parts of the correlation part 
of the self-energy for the minority spin state r' 25 . (b) The real and imaginary parts of the 
correlation part of the self-energy for the majority spin state [8] . 



Figure 7: NiO bandstructure. (a) Comparison between the LDA (solid line) and the experimental 
bandstructure (Shen et al 1990, 1991a,b). (b) Comparison between the LDA (solid line) and the 
quasiparticle bandstructure in the GWA [11]. 



Figure 8: The energy-loss spectra of NiO. The smooth solid curve corresponds to the calculated 
spectrum with virtually no gap in the LDA Hamiltonian and the dashed one with «5 eV gap. 
The other curve is the experiment [12] 



r ^45@ mt^M^n (2000^Jt)j (WD 

Figure 9: The spectral function of NiO at the T point for a Ni 3d state. The solid line corresponds 
to the case where the starting LDA Hamiltonian has virtually no gap and the dashed line to the 
case with «5 eV gap. (b) The same as in (a) but magnified 10 times [12] 



Figure 10: Diagrammatic expansion for the Green function to second order in the GWA and the 
cumulant expansion respectively. The solid line represents the noninteracting Green function 
Go and the wiggly line represents the screened interaction W. 



Figure 11: (a) The experimental spectral function for Na (dots). The solid line is a synthetic 
spectrum obtained by convoluting the density of states from a bandstructure calculation and 
the the experimental core level spectrum. BG is the estimated background contribution. The 
data are taken from Steiner, Hochst, and Hufher (1979). (b) The total spectral function of Na 
for the occupied states. The solid and dashed line correspond to the cumulant expansion and 
GWA respectively [13]. 



Figure 12: Feynman diagrams for the T-matrix (square): The wiggly and the solid line with 
arrow represent the screened interaction U and the Green function G respectively. 



Figure 13: Ni spectral functions at the X point for the second band [97]. 



Figure 14: The quasiparticle dispersions for r s — 4 corresponding to full self-consistency (solid 
line), partial self- consistency (dashed line), first iteration (dotted line) and the free electron 
(dashed-dotted line) [54]. 



Figure 15: The partially self-consistent spectral function A(k = Q.hkp, w) compared to that of 
the first iteration for r s — 4 [106]. 



Figure 16: The fully self-consistent spectral function A(k = 0.5/u/r,t<;) compared to that of the 
first iteration for r s = 4 [53]. 



Figure 17: The spectral function of the partially self- consistent self-energy, T = |ImS|/7r, for 
k — kp compared to that of the first iteration [106]. 



Figure 18: The real part of the partially self-consistent self-energy for k = kp compared to that 
of the first iteration [106]. 
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